Correlated Dirac particles and Superconductivity on the Honeycomb Lattice 
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We investigate the properties of the nearest-neighbor singlet pairing and the emergence of d- 
wave superconductivity in the doped honeycomb lattice considering the limit of large interactions 
and the t — J\ — J 2 model. First, by applying a renormalized mean- field procedure as well as 
slave-boson theories which account for the proximity to the Mott insulating state, we confirm the 
emergence of d-wave superconductivity in agreement with earlier works. We show that a small but 
finite J2 spin coupling between next-nearest neighbors stabilizes d-wave symmetry compared to the 
extended s-wave scenario. At small hole doping, to minimize energy and to gap the whole Fermi 
surface or all the Dirac points, the superconducting ground state is characterized by a d + id singlet 
pairing assigned to one valley and a d — id singlet pairing to the other, which then preserves time- 
reversal symmetry. The slightly doped situation is distinct from the heavily doped case (around 
3/8 and 5/8 filling) supporting a pure chiral d + id symmetry and breaking time-reversal symmetry. 
Then, we apply the functional Re normalization Group and we study in more detail the competition 
between antiferromagnetism and superconductivity in the vicinity of half-filling. We discuss possible 
applications to strongly-correlated compounds with Copper hexagonal planes such as In3Cu2VOg. 
Our findings are also relevant to the understanding of exotic superfluidity with cold atoms. 



I. INTRODUCTION 

Recently, graphene systems have attracted a consid- 
erable attention of experimentalists as well as theorists 
[TJ [2] . Graphene which consists of a single layer of carbon 
atoms forming a honeycomb lattice allows to realize in a 
condensed-matter system the Dirac equation, where elec- 
trons behave as massless Dirac fermions. The observation 
of massless Dirac fermions in monolayer graphene has en- 
gendered a new era of science and electrons in graphene 
embody a typical example of weakly interacting relativis- 
ts quantum systems [3][4] . The chemical potential can be 
tuned and hence it is possible to change the concentration 
of carriers, holes or electrons, opening the door for car- 
bon based electronics. It is also relevant to mention that 
artificial graphene has also been realized in cold atom sys- 
tems [5] , with photons [BHH] and using scanning tunneling 
microscopy techniques [9 . Doped graphene exhibits a fi- 
nite density of states which favors antiferromagnetic spin 
fluctuations [10] and then may lead to unconventional 
superconductivity. Experimentally, superconductivity in 
graphene has been induced by proximity effect through 
contact with superconducting electrodes [11] . This shows 
that Cooper pairs can propagate coherently in graphene. 

Several theoretical attempts have been made to de- 
scribe the emergence of superconductivity in graphene 
[T2HT6] as well as the formation of zero-energy states in 
the cores of vortices or at the boundaries [T7HT9] . In gen- 
eral, if these bound states appear in even number then 
they are not topologically protected and, for example, for 
small coherence lengths [18] their energies approach the 
Caroli, De Gennes, Matricon energy [20]. Uchoa et al. 



[T2] suggested that an extended s-wave SC phase may 
be realized at the mean-field level due to the peculiar 
structure of the honeycomb lattice. On the other hand, 
for a purely on-site repulsive Hubbard interaction U, as 
shown through the functional Renormalization Group 
(fRG) [13] , the nearest-neighbor spin exchange interac- 
tion J can lead to a d+id (d xy +id x 2_ y 2) superconducting 
state as a reminiscence of the superconductivity on the 
triangular lattice [21 . A similar fRG scheme has been 
applied on the square lattice [22-25 . The d + id super- 
conducting state has also been found using a mean-field 
theory on a toy model with singlet pairing between dif- 
ferent sublattices [14] [15] . A similar result has also been 
confirmed via numerical results based on a recently de- 
veloped variational method, the Grassmann tensor prod- 
uct state approach [26]. These theoretical investigations 
concern the situation close to half-filling. 

On the other hand, experimental techniques in doping 
methods [27[ [28] have allowed to approach the van Hove 
singularities, which corresponds to dope the graphene 
close to the M point of the Brillouin zone, i.e., for 3/8 or 
5/8 electron filling (which corresponds to doping ±1/8 
from the Dirac points; pristine graphene corresponds to 
1/2 filling). The logarithmically divergent density of 
states at the van Hove singularities (van Hove filling) 
unambiguously favors the appearance of d + id super- 
conductivity for weak repulsive on-site interactions, as 
shown from a perturbative RG approach [16] and a fRG 
framework [29j [30] . Rather unique on the honeycomb lat- 
tice is the degeneracy of the two d-wave pairing channels 

[H[3T]. 

In graphene, the Hubbard interaction is approximately 
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half the bandwidth which places this material in the 
intermediate-coupling regime. Below, motivated by the 
recent realization of strongly-correlated honeycomb lat- 
tice materials such as InsC^VOg [32-34 , we are rather 
interested in the stronger interaction regime which al- 
lows to realize Mott physics in the half-filled situation. 
A similar situation could be eventually reached using cold 
atomic systems on the honeycomb lattice [5 . The hon- 
eycomb lattice, which is a bipartite lattice, allows for a 
spin density wave order [35-38 . In L23CU2VO9, singly 
occupied 3z 2 — r 2 electrons of coppers contribute to an 
antiferromagnetic moment with S = 1/2. At quite inter- 
mediate values of the Hubbard interaction, using Quan- 
tum Monte Carlo (QMC) simulations and an accurate 
finite-size scaling up to 648 sites, Meng et al. [39] have 
recently reported the emergence of a spin liquid ground 
state in the range 3.5 < U/t < 4.3, which is characterized 
by a single-particle gap where band theory would predict 
a metallic behavior; see also Refis. 40 and 11 . The fin- 
gerprints of such a Mott phase without long-range Neel 
ordering have also been reported using cluster methods 
[42] [43] when increasing the size of the cluster unit cell 
[42] and in anisotropic lattices [44]. Sorella et al. [45] 
extended the QMC calculations up to 2592 sites and did 
not find evidence for this spin liquid phase region. On 
the other hand, the existence of a spin liquid phase, with 
a spin gap and Z2 symmetry, has been corroborated in 
the strong-coupling J\ — J 2 effective spin model on the 
honeycomb and square lattices [46 -48 . For large values 
of J r 2 , one may also expect a dimerized symmetry-broken 
phase [46] [49]. On the honeycomb lattice, in particular, 
this invalidates the possibility of an algebraic spin liquid 
with U(l) or SU(2) gauge theories at relatively moderate 
interactions [50, 51 . On the other hand, stable algebraic 
spin liquids on the honeycomb lattice with Z2 symmetry 
have been predicted for strong interactions [52H54] . The 
undoped compound L13CU2VO9 seems to yield a mag- 
netically ordered ground state [33[ [34] . Finding a spin 
liquid in two or three dimensions represents a consider- 
able challenge in condensed-matter physics [55H65] . It is 
also relevant to note that a spin liquid ground state has 
also been reported in two-dimensional Kagome [66] [67] 
and organic triangular materials [68-7T] in relation with 
theoretical developments [50] [52] I60H62] [65] [72H82] . 

In this paper, we seek to start from the Mott insu- 
lating and Neel ordered phase on the honeycomb lattice 
and dope the system away from half-filling with a few 
holes. Our main goal is to investigate the emergence 
of pairing and superconductivity within the framework 
of the t — J\ — J2 model, when including a finite (but 
small) next nearest-neighbor spin exchange interaction. 
A strong correlation view of the Hubbard model, through 
the t — J model [83] , was advanced by Anderson [56 , who 
conjectured the relevance of a spin liquid phase or Res- 
onating Valence Bond (RVB) phase as a result of the 
motion of the holes, destroying the antiferromagnetic or- 
der. The RVB state corresponds to a spin-gapped sin- 
glet state with no symmetry breaking. The doped spin- 



1/2 honeycomb lattice compound L13CU2VO9 might be 
a good candidate for the realization of such a physics 
through the t — J\ — J2 model [84 . The condensation of 
the holes (bosons) at low temperatures should result in 
a superconducting ground state. 

On the square lattice, following the Gutzwiller pro- 
jector point of view [85] [86] . this scenario has been 
pushed forward through a projected mean-field theory 
(the renormalized mean-field theory or Gutzwiller RVB 
theory) removing all components of the wavefunction 
with doubly occupied sites [87-93 , "slave-particle" ap- 
proaches [54 , 94 -107 and powerful numerical approaches 
[108141 18] . We shall also mention some theoretical pro- 
gess accomplished close to the Mott state [9TJ I119[ I121J - 
H25] . 

By applying the renormalized mean-field theory or 
Gutzwiller RVB theory [87-92 and "slave-particle" ap- 
proaches [54 , 94-107 , first we will show that on the hon- 
eycomb lattice and close to half-filling, the d ± id pair- 
ing order parameter favoring spin singlet between near- 
est neighbors is stabilized by the small J2 antiferromag- 
netic spin exchange whereas the extended s-wave pairing 
strength diminishes. For strong interactions, as noticed 
in Ref. 1126) the ground state taking d + id in one valley 
and d — id in the other does not break time-reversal sym- 
metry in contrast to the heavily-doped situation at 3/8 
or 5/8 electron filling at weak interactions [16] [29 , 30J and 
minimizes the total energy since the whole Fermi surface 
becomes gapped. The edge state from the d + id pairing 
in one valley is canceled by that from the d — id pairing 
in the other valley. This assignment turns out to be es- 
sential because the d + id order parameter in one valley 
vanishes in the other valley, allowing the Dirac spectrum 
to be gapless [126 . Using the fRG approach, then we al- 
low for the presence of antiferromagnetism at half-filling 
and study more rigorously the competition between su- 
perconductivity and antiferromagnetism in the presence 
of the J 2 term, following the scheme of Ref. [13] 

The remainder of the paper is organized as follows. 
In Sec. II, we introduce the model Hamiltonian, discuss 
the dominant pairing symmetries using the renormalized 
mean-field theory (RMFT) and the effect of a finite next- 
nearest neighbor spin exchange J2. We also comment on 
the possibility of stable Z2 (gapped) spin liquids at half- 
filling for not too small values of J2. In Sec. Ill, we 
present the theoretical framework, the main equations 
and the results. In Sec. IV, by applying fRG, we address 
the competition between antiferromagnetism and super- 
conductivity as a function of J2 and doping. In Appendix 
A, we compare our results obtained from the RMFT with 
those obtained within the U(l) slave-boson theory. 



II. MODEL HAMILTONIAN 

To capture the effect of strong interactions (or Mott 
physics at half-filling) in honeycomb lattice compounds 
such as In3Cu2VOg [32-34 , we consider the (renormal- 
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ized) t — J\ — J2 model: 



H = -tg t (4* d j° + h - c 



(1) 



+ Ji9i ^2 Si • Sj + ^2^2 Si ' 

(iJ) ((i,3)) 

The Gutzwiller projector [85 ensuring that configura- 
tions with doubly occupied sites are forbidden is replaced 
by statistical weighting factors, g t — 25/(1 + 5) |127j and 
9i — 92 = 4/(1 + 5) 2 [87 , and explicitly depend on the 
doping level 5 [57H92] . Note that here, 5 = 1 — n where 
n refers to the number of electrons per site. For the half- 
filled situation, the number of holes per site is 5 = 0. 
In addition, denotes a nearest neighbor pair where 
i < j is assumed. We also denote c and d electron anni- 
hilation operators associated with the two sublattices (A 
and B) of the honeycomb lattice. 

In the limit of large on-site interaction, this model can 
be derived from the Hubbard model, similarly to the 
derivation of the Kondo model from the Anderson model 
|128] , by resorting to a perturbation theory mt/U up to 
fourth order processes [46j [~ 



Ji 



4*2 
U 
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J2 



4^ 
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Note that an antiferromagnetic phase has been reported 
for J2/J1 < 0.08 [46 which corresponds to U/t ~ 4.3. 
This is in agreement with QMC results of Refs. [39] and 
l40l and Cluster methods of Refs. |42| and |43[ but in con- 
trast with the very recent QMC results found in Ref. 
l45l which predict that the spin density wave order on the 
honeycomb lattice would appear simultaneously with the 
Mott transition at lower interaction strength. 

Below, we start from half-filling with the spin density 
wave order where J2 < O.O8J1. Note that for the undoped 
compound L13CU2VO9, it has been recently estimated 
that J2/J1 ~ 0.04 [84 . The undoped compound seems to 
order at low temperatures [32H34] . Hereafter, we do not 
focus on the half-filled situation and assume that there 
is a finite hole concentration. By increasing the number 
of carriers, one may expect an RVB type scenario and a 
gapped spin liquid [56], as a result of the motion of the 
carriers (holes). Hereafter, we describe this aspect of the 
problem through the t — J\ — J2 model close to half-filling. 



A. Pairing symmetries 

Firstly, we investigate the t-J model where J\ = J and 
J 2 = applying the RMFT |87U92l [T30] . The emergence 
of nearest neighbor singlet pairing for repulsive on-site 
interactions on the honeycomb lattice has been first dis- 
cussed in Refs. [13] and [14] Our procedure is slightly 
different from the one used by Black- Schaffer and Do- 
niach [14] since we take into account the large interaction 



limit through the statistical weighting factors g t , g\ and 
<72> an d we shall also address the effect of a finite next 
nearest neighbor coupling J2. 

Following the procedure used on the square lattice [87- 
[92], it is convenient to introduce the mean-field order 
parameters, 



Xij = ^9lJ ^(ctrdja) 



(3) 



and we focus on the nearest- neighbor singlet pairing con- 
tribution (on-site pairing is forbidden due to the very 
large on-site repulsion). Assuming the uniform solution 
for the x field the mean-field Hamiltonian takes the form, 



H = (-tg t - f ) E ( c Uv + 4 C -) ( 4 ) 

-\T.{^{ c lAi+ d )Ai)+ h - c ] 

(ij) 
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Here, we have released the constraint i < j and intro- 
duced the chemical potential \i explicitly. It is judicious 
to Fourier transform the Hamiltonian and introduce the 
symmetric and antisymmetric (band) combinations of 
the electron operators c and d which diagonalize the ki- 
netic part, as follows 



d \ta = ^= exp(-Z0 k )(/ ko - - #ka). 

This results in the Hamiltonian: 

ka kcr 
k 

+ E( A k (flW-U-9tA-u)+ h - c -) 

k 



(5) 
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Jgi 



(6) 



Hereafter, the sum over a corresponds to a summation 
over the three nearest neighbors on the honeycomb lattice 
and A a is defined in Eq. (7). Here, N s corresponds to the 
total number of sites, is the intraband pairing while 
A£ is the interband counterpart breaking time reversal 
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symmetry, 



At 



a 

Q! 

ek = (-tfft - |) |7k|- 



0k) 



(7) 



ik-R Q 



The phase 0k 
k and satisfies 



It is convenient to define 7k = ^2 a e 
is defined as 0k = arg(^ a e zk ' Ra ) 
the following relation exp(z0k)7k = exp( — £0k)7k = l7k|- 
At a general level, the intraband pairing contribution 
exhibits an order parameter even in k space, and one can 
check that it corresponds to the singlet pairing state. In 
contrast, the interband pairing contribution has an or- 
der parameter odd in k space. For a bond-independent 
s-wave order parameter, the interband (spinon) pairing 
then is identically zero, but this is not necessarily the 
case for an arbitrary wave symmetry (such as d-wave 
symmetry). In this paper, we mostly consider the two 
dominant nearest-neighbor singlet pairing order param- 
eters when assuming purely on-site interactions, namely 
the extended s-wave pairing (ES) and the d±id pairing 
p~3| Hi] . For the ES pairing only the intraband pairing 



form factor are non-zero, 



l£ Ael k.K Q 



(8) 



For the d±id pairing state, the order parameter can be 
viewed as a mixture of d xy and d x 2_ y 2 [T3| [14] : 



,d±id 



(|) A x 2_ y 2(k) ±isin (|) A^(k). (9) 



In fact, it is perhaps judicious to remember that the 
d x 2_ y 2 and d xy wave symmetry functions satisfy: 



d>x 2 — y 2 {k x i ky) — ^ 



J V3 



d X y(k x ,k y ) = ie sin 



e ^ cos 



(10) 
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Therefore, it is convenient to introduce the notations: 
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V 1 

1 \r 



^ ^ A a cos(k • K a - k ) = 

ex. 

l^ e Ma-D/3 cos(k . Ra 

a 

i ^2 A ai sin(k • K a - k ) 

a 

I^ e 2M«-D/3 isin(k . Ra 



ATI 
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FIG. 1. Form factors of the ES and d + id pairing solutions. 
Due to the overlap between the Fermi surface and the nodes 
of the ES solution, the ground state should favor the d + id 
pairing as discussed in Sec. III. On the other hand, to gap 
all the Fermi surface close to half-filling, the solution that 
minimizes the whole energy will be taking the d + id pairing 
solution in one valley and the d — id in the other [126 . 



could anticipate that the ES solution will have a higher 
free energy then favoring the d ± id pairing symmetry. 
The mean-field equations will be discussed in the next 
Sec. III. 



B. Effect of J 2 

To investigate the effect of the next nearest neighbor 
spin exchange J2 on the physical properties of the sys- 
tem, here we will assume that since we consider the limit 
where J\ ^> J 2 , the dominant pairing order parameter is 
Aij which corresponds to (d-wave symmetry for) nearest 
neighbor singlet pairing. As a result, we only introduce 
an extra particle-hole order parameter which couples the 
next nearest neighbor sites: 



X 



(12) 



The mean-field Hamiltonian on the honeycomb lattice 
then takes the following form: 



H 



= (~tgt - 



(13) 



X_ 
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E 



(ij) 

1 v Wi 

3 frrf Jm 



I A;, | 2 



3 j-i Jm 



h h.c. 

N s z' |xf 
3 J292 



( C la C ^ + 4v d i<r 



In fact, owing to the large overlap between the nodes of 
the ES form factor and the Fermi surface (see Fig. 1), one 



Here, z' denotes the next nearest-neighbor coordination 
number, i.e., z r = 4 for the square lattice and z r = 6 for 
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the honeycomb lattice. After Fourier transformation, we 
obtain: 

H = (-tgt - |) ^bkcLA, + 7k4,Ck4 (14) 

kir 

x' \ 

L ka 

- \ E{ A k(4/- k ; - c U d Ut) + h.c.) 

k 

| ly\ X ? , ly |A^| 2 | 7V^|xf 
3 ^ Ji#i 3 ^ Ji^i 3 J 2 g 2 



kcr 



where, 



7k = E eikRa 
Ck = a = E e ' k ' R ' 3 



(15) 



(16) 



The sum over f3 here denotes the summation over next 
nearest neighbors. 

Assuming that J\ ^> J 2l we then have a modified dis- 
persion relation for the fermions £k = (—tgt — f)l7 k l — 
^-Ck and an additional constant term 2N ^ X I . The self- 
consistent equations for the t — J\ — J 2 model then will be 
directly inferred from the ones for the t — J model. How- 
ever, we shall have an additional equation determining x' • 
In Sec. Ill, we shall show that the J 2 term through the 
extra order parameter x' wm hdp stabilizing the d±id 
spin pairing. 



III. RESULTS FROM RMFT 



Hamiltonian can be formally re- written as: 

H = E + Y1 E u{ a li a u + a -u a -u} ( 17 ) 

J191 J\g\ 

Here, the sum I = 0, 1 stems from the path integral of 
the two-band Hamiltonian and 



hi 



kZ 



^ ki = \/{(-i)'a-M} 2 + i iA k | 2 , 



(18) 



where we have introduced E^i = y + ||Ak| 2 with 

A k = E tt Ae M « and fa = (-l) l (-tg t - §)| 7k | - M 
when J 2 = 0. The free energy then takes the form 



F = -2T^ln (2 cosh 
kz ^ 



k/ 



(19) 



For simplicity, we set the Boltzmann constant ks — 1- 
At the stationary point of the free energy F we obtain 
the BCS-like self-consistent equations, 



tanh 



kZ 



N s ^ E kl 
X = -^V ( - 1) >' l7k 'tanh-^ 



(20) 



i.y„ ^ £ ki 
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The solution of these equations will be discussed in Sec. 
III. C. Now, we present the equations for the d±id case. 



Here, we give the main equations associated with the 
different pairing solutions by resorting to the renormal- 
ized mean-field theory. 



Extended s-wave scenario 



B. d-wave scenario 



For the d±id case, we can rewrite the Hamiltonian by 



introducing $1 = [f^, /_ k; , g^, #_ k J, 



As discussed earlier in Sec. II, for the ES pairing, only 
the intraband pairing form factors of Eq. (9) are non- 
zero. Since the effect of J 2 is relatively simple follow- 
ing the scheme above, we present the main equations for 
J 2 = 0. After standard Bogoliubov transformation, the 

I 

Thus we can determine the energy dispersion of the Bogoliubov quasiparticles, 
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FIG. 2. Spin gap Ad+id and A s at T = (denned in Eqs. (8) and (11)) and Superconducting Transition Temperature 
(T c ~ gtA) when J2 = for the d±id and ES scenarios as a function of the hole doping parameter S = 1 — n; the half- filled 
case here corresponds to one electron per site or (5 = 0. The order parameters are taken in units of 3g s J/4. Inset: one can 
observe that the particle-hole order parameters behave identically for the d±id and ES situations. 



and for the d±id case we use the ansatz in Eq. (12). The self-consistent equations then read, 

dE kl = A\r k \* | A|ij| 2 | 1 (-i)'(i6A 3 |rLI 2 |r^| 2 + 2A(ao-Cki) 2 |r^| 2 ) 

OA E kl E kl 4E kl ^/T^T _ £ i)2 + A 2 |lt| 2 (ao - 6a) 2 + 4A4|rj c | 2 |r£| 2 
d_E^ = J^ + J: 

dfi Ekl 2Ekl ^JMo - &) 2 + A 2 |r£| 2 (£ k0 - ai) 2 + 4A4|p fe | 2 |ri| 2 (21) 

ggkj = -(gko-^ki)lTkl , 0.5(-l)'+ 1 (^ -^ 1 )(eko + eki)|7k|-(-l)'4|A| 2 |r^| 2 | 7k |e k 

dx AEu 4i? ki - CL) 2 + A 2 |r^| 2 (a - &i) 2 + 4A4|ri c | 2 |r^| 2 ' 

Note that near the superconducting transition temperature this results in 

dEu A\Tl\ 2 A(ao + &i) + (-l)'A(C ko -&i) A|4| 2 _A . 2 

3A £ ki + 3u(&o + &i) ' k ' E kl - M ' k ' ' 



Finally, the self-consistent mean field equations for the 
d±id situation take the form 



S = — — t~ tanh 



A = ^y t anh- 

27V,- ^ 



2iV. 



rim tanh 



k; k; 



2 <9A 

2 dx ■ 



(22) 



Taking the limit = is also consistent with the 
equations for the ES case discussed above. 



C. Spin Gap and Superconductivity 

As a result of the strong on-site repulsion, the domi- 
nant pairing term is between nearest neighbor sites and 
counting the lattice symmetry then the natural candidate 
will be d x 2_ y 2 or d xy for nearest neighbor singlet pairing 
as a reminiscence of the triangular lattice [21] . At a gen- 
eral level, one can introduce a general combination of 
d x 2_ y 2 and d xy wave pairing for the pairing term: 

A k = cos0A x 2_ y 2(k) ±ismOA xy (k). (23) 

As already elaborated in Ref. [14] one can show that 
the minimum of free energy occurs for — tt/3. Such 
a pairing solution could favor a stable Zi spin gapped 
phase at half- filling for not too small J2, as supported 
in Refs. [46] [48j 126 . The Z2 symmetry can be seen 
from the mean-field Eq. Q following similar arguments 
as in Refs. [48] and 11261 It is perhaps important to un- 
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FIG. 3. Spin gaps (at T = 0) versus doping, for J 2 / 0, from 
the RMFT. Notations are similar to those in Fig. 2. 



derline that (gapped) spin liquids are stable in two 
dimensions beyond mean-field arguments [131] in con- 
trast to certain U(l) analogues [38| |5Q| [T32] . Other ver- 
sions of spin liquids might be protected from gauge fields 
in the large N limit [133-136J. On the other hand, for 
J 2 — >• 0, the ground state at half- filling is an antiferro- 
magnet [38j [40j [45] , and the proximity to antiferromag- 
netism will be addressed thoroughly through the fRG. 

In Fig. 2, we present our results for the pairing 
strengths of the ES and d ± id situations close to half- 
filling when J 2 = 0. Within the RMFT, the quantity g t A 
can be interpreted as an "approximate" superconduct- 
ing transition temperature [87-90 (remember that this 
approach ignores the possibility of antiferromagnetic or- 
der). As already anticipated earlier, we confirm that the 
d ± id solution is more favorable than the ES scenario. 
For completeness, we compare our results for the spin gap 
(RVB gap) and superconducting transition temperature 
at J 2 = with those obtained within the U(l) slave bo- 
son approach adapted to the honeycomb lattice; consult 
Appendix A for a comparison with the slave-boson the- 
ory. Results obtained via the slave-boson theory are in 
qualitative agreement with the RMFT. By doping with 
holes, in the strong coupling limit, we then confirm the 
occurrence of a d ± id superconducting ground state; on 
the other hand, a more refined (probably numerical) ap- 
proach would be necessary to estimate the evolution of 
the ground state to the heavily doped case in the case of 
strong interactions. Let us emphasize that for weak in- 
teractions, a d±id superconducting ground state break- 
ing time-reversal symmetry has been found close to 3/8 
filling (S = 1/4) [EH [29]. 

As shown in Fig. 3, the relative strength of the ES 
gap becomes less pronounced when including the effect of 
the finite next nearest-neighbor spin coupling J2, making 
the superconducting transition towards the d±id ground 
state more favorable. On the other hand, the RMFT ig- 
nores the presence of long-range antiferromagnetism and 
therefore the competition between superconductivity and 
antiferromagnetism will be studied via the fRG. 



IV. RESULTS FROM FRG 

Here, we want to describe what information can be 
gained beyond the RMFT by using a fRG analysis of 
the unconstrained J1-J2 model. To this end, we adapt 
the fRG approach for interacting fermions in the the so- 
called TV-patch approximation (for a recent review, see 
Ref. 137) to study the leading instabilities of the J1-J2 
model on the honeycomb lattice; see Fig. [4j 

A. Methodology 

The fRG treatment offers a) an unbiased comparison of 
the different possible instabilities or ordering tendencies, 
b) provides estimates for energy scales of these instabili- 
ties, both including the coupling of different fluctuations 
beyond mean-field theory. Note however, in contrast 
with the previous sections, we study the unconstrained 
model and do not use the Gutzwiller projection. The 
reason for this difference is that the fRG is a technique 
that is perturbative in the interactions, and therefore the 
weakly doped situation in the Gutzwiller approach with 
the small renormalized hopping term ~ S is not a good 
starting point for this method. So, in principle, the fRG 
approach for the unconstrained J\ — J2 model does know 
about the antiferromagnetic spin interactions on neigh- 
bored sites, but not about the strong onsite correlations. 
In order to make up for this, we also include a moderate 
local Hubbard interaction and check whether our results 
depend qualitatively on this. 

The fRG scheme employed here is the same was as 
recently used to explore mono- Q3J [138], bi- [29l fT39l fT4Q] 
and trilaver [l41] honeycomb models with density-density 
interaction terms. A Brillouin zone discretization in N 
angular patches around the K and K' points is em- 
ployed in order to resolve the wavevector- and band- 
dependence dependence of the scale-dependent interac- 
tions Va(A?i, ni, ^4) (where the ks denote 
the incoming and outgoing wave vectors, and the ns the 
band indices of the interaction). Upon integrating out 
the electronic degrees of freedom with lowering an in- 
frared cutoff energy scale A, the wavevector- dependent 
one-loop corrections (particle-hole and particle-particle 
bubbles) to the bare interactions are summed up to infi- 
nite order. The standard approximation employed in this 
instability analysis are, just as in many previous studies 
(e.g. cited in Ref. 1137)) . that the electronic self energy is 
ignored, and that vertices of order higher than four are 
not taken into account. Recent work on the square lat- 
tice Hubbard model has shown that the inclusion of the 
self energy does not change the conclusions |142[ 1143] . 
For sufficiently strong interactions, this leads to a flow 
to strong coupling at a nonzero cutoff scale A c where 
a part of V/y(fci, ni, &2> &3> ^3? n*) seems to diverge. 
The scale A c can be used as an estimate for the energy 
scale of the ordering phenomenon suggested by this flow 
to strong coupling. Furthermore, the wavevector- and 
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FIG. 4. Left panel: Interaction vertex labeled with the 
spin convention (upper diagram). Below, the loop contri- 
butions to the flow of the interaction vertex including the 
particle-particle diagram (a), the crossed particle-hole dia- 
gram (b), and the direct particle- hole diagram (c). Right 
panel: U N = 24" -patching scheme of the Brillouin zone with 
constant wave vector dependence within one patch and the 
representative wavevector chosen on the Fermi-line shown 
here for three different choices of the doping S. The red dashed 
line corresponds to van Hove filling. 




0.00 ' 0.05 0.10 0.15 0.20 0.25 
S 



FIG. 5. fRG critical scale A c in units of the hopping t vs. 
doping 6 at T = for various choices of parameters J\ , J2 
and onsite interaction U. The point 5 = 0.25 corresponds to 
van Hove doping. J2 > suppresses the AF SDW ordering 
tendencies, but has only smaller quantitative effects on the 
(i-wave pairing at larger couplings. Inclusion of U does not 
change the qualitative findings. 



band-dependence of the leading terms in this divergence 
allows one to extract the order parameters that may ac- 
tually order below this instability scale. 



B. Results 

First let us study the case with pure spin-spin interac- 
tions and set U = 0. We explore a region of the phase dia- 
gram near the charge neutrality point and with J\ ^> J 2 . 
In Fig. [5] we plot the evolution of the fRG critical scale 
A c as function of the density deviation from half filling 
for Ji = 1.6t and three choices of small J2. In the plot 
we also indicate the leading instability (for a description 
how these phases are identified from the running cou- 
plings, see, e.g., Ref. 13). For zero and small doping, 
we get an antiferromagnetic (AF) spin-density-wave in- 
stability (SDW) at rather large scales. Note that due 
to the Dirac cones in the dispersion, a nonzero minimal 
interaction strength is required to obtain an instability 
at half filling. The critical value comes out as J\ ~ lAt. 
From the experience with the honeycomb Hubbard model 
with pure onsite interactions (discussed e.g. in the bi- 
layer case in Ref. I140p , we expect that the fRG in this 
approximation will somewhat underestimate the minimal 
true value, but qualitatively get the correct picture. The 
AF-SDW state was also found in QMC studies of the 
larger-/7 honeycomb Hubbard model [40l [4TJ |45] . Hence, 
our study without onsite correlations fits in quite consis- 
tently. 

Regarding the effect of J 2 > 0, in Fig. [5] it can be seen 
clearly that J2 reduces the AF-SDW tendencies quite 
strongly. This frustration effect is of course not unex- 
pected, and is truthfully captured by the fRG that allows 
for a coupling of fluctuations with different wave vectors. 

When we increase the doping, the critical scale for the 



AF-SDW instability drops strongly. Then, beyond a crit- 
ical doping which depends on J2, a pairing instability in 
the d-wave channel takes over, at doping levels which 
can be inferred from Fig. |5j Here, for symmetry rea- 
sons, the pair scattering in the cL^-channel and in the 
d x 2 _ y 2 -channel diverge together. In Fig. [6] we show snap- 
shots of the wavevector-or patch-dependence of the effec- 
tive interactions near this instability, transformed back 
into the sublattice basis, for different combinations of 
the sublattice indices. In the left panel, incoming and 
outgoing particles are all on the same sublattice. There, 
no strong interactions can be observed, i.e. the effective 
interaction does not have any large intra-sublattice con- 
tributions. The picture changes in the middle and right 
panel, where strong diagonal features with large positive 
and negative interactions can be found. These sharp lines 
occur for incoming wave vectors (labelled by the patch 
indices k\ and A^) adding up to zero, i.e., they belong 
to the Cooper pair scattering channel, and the instabil- 
ity should be interpreted as Cooper pairing instability. 
The sign structure along this line encodes the symme- 
try of the Cooper pair (p, —p), when p moves around the 
Fermi surface. In order to see this symmetry more clearly, 
we plot in Fig. [7] the pair scattering k,—k — » p, — p 
with p varying around the Fermi surface in the Bril- 
louin zone hexagon, with k held fixed near the Brillouin 
zone boundary near K, all in the band which crosses the 
Fermi level. We clearly see the modulation of the pair 
scattering with p. We also plot the d-wave form fac- 
tor V d (k,p) = -Vo d% y (k)d xy (p) + d* x2 _ y2 (k)d x 2_ y 2(p) . 

We can see that the pair scattering in the effective inter- 
action near the instability follows this form factor rather 
well, both for J2 = and for nonzero, small J2. For 
comparison, we also plot the form factor for extended 
s-wave pairing on nearest neighbors. This does not give 



9 





50 



FIG. 6. Typical effective interaction vertex near the critical 
scale in the regime of the d-wave instability in units of t. Left 
Panel: Orbital combinations with 01=02=03= 04 where 
Oi G {a, b}. The numbers on the axis specify the number 
of the patch as shown in Fig. [4] On the horizontal axis the 
wavevector k\ can be read off and on the vertical axis we 
enumerate fe. k$ is fixed on the first patch, k± then follows 
from momentum conservation. Middle Panel: Effective vertex 
function for the orbital combination, where o\ — 03, 02 = 04 7^ 
01. Here, we can clearly identify sharp diagonal structures 
(k\ — — fe) with a d- wave- modulation of the amplitude along 
the diagonal, see Fig. [7] Right panel: Effective vertex function 
for the orbital combination, where o\ — 04,02 = 03 7^ o\. 
Also for this orbital combination a sharp diagonal structure 
emerges. 



any good match for the fRG data and confirms the strong 
dominance of the d-wave pairing tendencies. 

In the previous sections, based on the RMFT, it was 
argued that the most stable pairing state in presence of 
the two degenerate d-channels would be to switch from 
d + id at one Fermi pocket and d — id at the other. The 
energy benefit from this comes due to the full gap now 
open on both Fermi circles, while a stiff d + id or d — id 
through out the BZ would have gap minima on one of 
the circles. If the pair scattering between the two Fermi 
circles (e.g. k at K and p at K') is rather weak, the 
pair scattering will not cause a sufficient energy penalty 
to prevent this switching the phase of the d-wave super- 
position from one Fermi pocket to another. With the 
fRG approach, without extended subsequent mean-field 
study of the low-energy model like e.g. in Ref. I144| it is 
not possible to make any refined statements about what 
would be the best paired state. Note however that the 
pair scattering from the fRG very closely follows the sim- 
ple nearest-neighbor form factor that also shows up in the 
mean field theory. In particular, the inter-pocket scatter- 
ing between the two Fermi circles is small and of varying 
sign, while the scattering within one circle is stronger. 
In fact, in the fRG data the inter-pocket scattering is 
even weaker than for the nearest-neighbor form factor, 
which would enter the mean- field treatment. Hence, the 
prerequisites for switching the phase of the d + id linear 
combination from one Fermi circle to the other to d — id 
are all there, and the fRG supports this energy lowering. 

Note that here we do not go to larger doping where the 
Fermi circles get close to each other or open in the middle 
between the K and K' points. Here we expect that the 
near-decoupling between the K and K' point in the pair 
scattering is not valid any more, and that a unique linear 



-50 



50 



10 20 

around hexagon 



-50 



OOOq poop 

^ V 



0006 



10 20 

around hexagon 



FIG. 7. Pair scattering Va(/c, — k — >> p,—p) near the insta- 
bility for doping S = 0.15, with k fixed one one discretiza- 
tion point near the zone boundary, and p moving through the 
other points around the hexagon. The circles are the fRG 
data, the dashed line is the nearest-neighbor d-wave form fac- 
tor oc — [d* 2 _ y 2 (k)d x 2_ y 2 (p)-\-d* y (k)d xy (p)], and the solid line 
the extended s-wave form factor. The left plot is for J\ — 1.5t 
and J2 = 0, the right plot for J± = 1.5t and J2 = 0.01£. 



combination of the degenerate d-basis functions needs to 
be chosen, as shown in Refs. [T6lor [29l 

Let us now discuss the effect of J2 on the pairing insta- 
bility. As can be seen in Fig. |5j the critical scale A c for 
d-wave pairing at small S drops when J2 is raised from 
0. A natural explanation of this observation is that the 
d-wave pairing is AF-spin-fluctuation driven and hence 
the reduction of the AF-SDW tendencies by J2 > also 
reduces the pairing tendencies. However, the character 
of the leading instability remains unchanged by this scale 
change, i.e. is still of d pairing type. In fact, the pair scat- 
tering near the instability follows the nearest-neighbor d- 
wave form factor even more closely than for J 2 = 0, as 
can also be seen in Fig. [5] We think that this is due to the 
lower scale, which implies less competition between the 
remnants of the SDW tendencies. Again, the extended 
s-wave pairing channel does not appear to be relevant, in 
agreement with the RMFT. 

Furthermore, including a nonzero U > in order to ac- 
count for local correlations does not change the character 
of the instabilities drastically. This has to be expected, 
as also the local repulsion drives SDW tendencies, so it 
basically adds to the nearest-neighbor interactions J\. 
Correspondingly, a smaller J\ ~ O.St can be used to pro- 
duce similar critical scales for the SDW as for U = (see 
Fig. [5]). On the other hand, the d-wave pairing instabil- 
ity occurs at somewhat lower scales than for U = 0. This 
ties in with the earlier observation that a [/-interaction 
alone with J\ — does not lead to pairing instability at 
reasonable scales in the slightly doped case [13]. 



V. CONCLUSION 

From the fRG study of the J\ — J2 model on the hon- 
eycomb model, we state that the weakly coupled model 
exhibits a standard AF-SDW instability above a criti- 
cal J\ when the doping and the frustrating J2 are not 
too large. Doping further in this regime leads to well- 
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formed <i-wave pairing instability, with predominant pair 
scattering within the respective Fermi circles around K 
or K' . Although the £RG approach dos not employ any 
Gutzwiller type renormalization and local correlation ef- 
fects enter only perturbatively through the Hubbard U, 
these findings tie in very consistently with the RMFT 
approach presented in the other sections. Our results 
suggest that the slightly doped compound InsC^VOg 
[34] could reveal an RVB phase as well as a (high-Tc) 
superconducting phase where the ground state is charac- 
terized by a d + id singlet pairing assigned to one valley 
and a d — id singlet pairing to the other, which then pre- 
serves time-reversal symmetry. The RMFT predicts that 
the T c might be quite high not too close to half-filling, 
as illustrated in Fig. 2. This analysis, that takes into 
account the J2 term, could be extended to multi-layer 
systems [TS3 H39HT4T1 [T45] . We could also explore the 
effect of next-nearest neighbor interactions [13]. 
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Appendix A: Slave-Boson Theory 

Here, we provide details concerning the U(l) slave- 
boson theory which has been used to give Fig. 3. Ex- 
tending the usual slave-particle procedure [104] on the 
honeycomb lattice, we represent the (slave) fermionic 
(here, spinon) operators on sublattice (A,B) respectively 
via {<4 a , c^}, {d\ a ,dka} (to build a connection with 
the RMFT) and the related bosonic (charge) operators 
{aj, a.j}, {b\, bi}. The mean-field Hamiltonian reads: 

H = 5Z( ekC L^ka + h.c) ~ fl c ^( C L C ko- + 4<Aa) 
k,cr k,<7 

+ ^(^k^^k + h -°) - ^ y^frjtfrk + a k a k) 

k k 

- E( A k( C k/-ki " 4/-k t ) + h.C.) 
k 

+N s (Qt XbXc + ^ + ^(l-<5) 2 + ^E A l A ^ ~ 2A ) 



(ij) 



(Al) 



where again N s denotes the number of sites on the lattice 
and we have defined the order parameters (which have 
been chosen to be slightly different from those introduced 



within the RMFT): 

A h = (4A - 4 A 

Xc = Xd = ^2(4a C 3 



Xb = Xa = (b%), 



(A2) 
(A3) 

(A4) 



and defined (Rj — R^ = R a ) 



ek = 


J \ 

{-txb - 2 Xchk 




Q J 


CJ k = 


-tXdk 




-A 


A k = 





(A5) 



We have introduced a Lagrange multiplier A to reinforce 
the condition of half-filling, for example on sublattice A: 
\(a\di + c\ a Ci a — 1) and similarly on sublattice B. 
We can now proceed and diagonalize the problems by 
introducing a transformation for the spinons similar to 
Eq. (5) in the main text. An identical procedure is then 
applied to the bosons (chargons). 
The free energy takes the form: 

F = -- R Y. ln ( 2 cosh ^H) " \ E ln ( X - e«"> +Wkl) 



k,z 



kJ 



■Z)-^/ + (-i)^k 



kJ 



+N s {6t XbXc + + ^(1 - Sf + W AjjAy - 2A} 



(A6) 



I = {0, 1} stems from the path integral of the two-band 
Hamiltonian. By analogy with the RMFT, we set 



E kl = x /& + |A k p 



J 



6d = -He + (-l)'kkl = -Mc + (-mtXb + 2 Ac;i7k 
giving the self-consistent equations for the ES case: 



^^tanh^ 
W s j^ E kl 2 



a = — E 



-E 



(-l)'Cki|7k| 



k,Z 



tanh 



kl 



(A8) 



Xb 



(-iyi7ki 

s k,J 



— E 

1 AT. ^ 



The chemical potential for fermions fif and bosons A 
are decided by ^ (flf ia ) = 1 - 6 = - ^ §^ and §f = 
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FIG. 8. Phase diagram obtained with the U(l) slave-boson 
approach for the case J2 = (here, J/t = 1). Details of the 
theory are presented in Appendix A as well as the definitions 
of the temperature scales Tbe and Trvb which correspond 
respectively to the temperatures associated with the boson 
(chargon) condensation and spin gap formation, respectively. 
Within this approach, an upper bound on the Superconduct- 
ing Transition Temperature is given by Min (Tbe, Trvb)- 



(actually is also determined by /i, implicitly. The 
Lagrange multiplier A connects the S = (b\bi) and /i), 



1 



^ tanh ~ 



(A9) 



and 



2N* ^ e^( w w+Mb) 



(A10) 



(biWA^O, and (fyi-fl&)?0, 



To evaluate the superconducting order parameter 
assume (bibj) 

i. e., we can numerical solve the self-consistent equations 
to find the temperature Tbe f° r holon condensation and 
Trvb for spin gap formation; see Fig. 8. 

We have checked that the self-consistent equations in 
Eqs. (A8) are consistent with those obtained within the 
RMFT in Eqs. (21) in the main text. Furthermore, we 
straightforwardly obtain similar equations as Eqs. (22) 
(in the main text) for the d±id situation. 
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